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We study the phase-synchronization properties of systolic and diastolic arterial pressure in healthy 
subjects. We find that delays in the oscillatory components of the time series depend on the 
frequency bands that are considered, in particular we find a change of sign in the phase shift going 
• from the Very Low Frequency band to the High Frequency band. This behavior should reflect a 

' collective behavior of a system of nonlinear interacting elementary oscillators. We prove that some 

models describing such systems, e.g. the Winfree and the Kuramoto models offer a clue to this 
phenomenon. For these theoretical models there is a linear relationship between phase shifts and 
' the difference of natural frequencies of oscillators and a change of sign in the phase shift naturally 

emerges. 

PACS numbers: 05.10.-a,05.45.-a,87.19.Uv 



I. INTRODUCTION 



Time series of physiological origin very often display synchronous behavior. Most likely, this is the result of collective 
behavior of a huge number of nonlinearly interacting elementary oscillators. Different examples of this phenomenon, 
as well as models of it, can be found for example in In the present paper we address a further example, i.e. 
I the relation between diastoHc (DAP) and systoHc (SAP) blood pressure signals in healthy subjects. The analysis 
' of blood pressure fluctuations may provide significant information on the physiology and pathophysiology of the 
Tj" ' autonomic control of the cardiovascular function The synchronization of these signals is expected, though a 

detailed study of its features is apparently still lacking. In a previous paper jjj it was noticed that DAP and SAP are 
ty^ ' characterized by a phase lag in the very low frequency band (VLF). The analysis of Q uses Fourier analysis, which is 
O , not particularly useful when non-stationary effects play a relevant role. In the present work we address two questions 
c/3 ' about the DAP/SAP relationship: Is the phase lag depending on the frequency band? Is the phase lag connected to 
^~>i a causal relation between SAP and DAP? To address these questions, we measured DAP and SAP signal in a number 
of healthy subject. Studying the mutual interdependency between the two time series, we conclude that there is not 
a causal relationship between DAP and SAP time series, i.e. none of the two is driver for the other. Moreover, a 
significant phase delay is found, for healthy subjects, in the VLF band and in the high frequency (HF) band. The 
phase shift between DAP and SAP is positive in VLF band and negative in the HF band. This change of sign in the 
phase shift has its origin in the regulatory mechanisms of blood circulation. A physiological interpretation of these 
■ mechanisms is beyond the scope of our work; however the hypothesis that synchronization results from the collective 
behavior of elementary nonlinear oscillators may offer a clue to its understanding. To exploit this idea we use below 
two well known models of coupled oscillators, the Winfree model and the Kuramoto model 0]. 

Winfree's paper jj] on coupled oscillators provided one of the first tractable examples of a self-organizing system. 
Winfree introduced an approximation that has become the standard approach to the study of ensembles of biological 
oscillators: In the weak coupling limit, amplitude variations could be neglected and the oscillators could be described 
only by their phases along their limit cycles. Winfree also discovered that systems of oscillators with randomly dis- 
tributed frequencies exhibit a remarkable cooperative phenomenon, reminiscent of a thermodynamic phase transition, 
as the variance of the frequencies is reduced. The oscillators remain incoherent, each running near its natural fre- 
quency, until a certain threshold is crossed. Then the oscillators begin to synchronize spontaneously. Winfree model 
was subsequently modified by Kuramoto who provided an analytically solvable version of it 0| . This field of study has 
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been very active all along and the analysis of synchronization phenomena remains a thriving area of research, see [7| 
for a review. Having in mind our experimental findings on the SAP/DAP relation, we first examine in Section UTI the 
phase shift between coupled oscillators in these models, once synchronization has been reached. We observe that there 
exists a simple linear dependence between phase shifts of synchronized oscillators and the difference between their 
natural frequencies. This phenomenon, to our knowledge never noticed before, offers a simple mechanism to describe 
the change of sign in the phase lag as the frequency band is changed. In Section IIIII we describe the experimental 
data and analyze them using the theoretical approach of Section Moreover we discuss the problem of the causal 
relation between the DAP/SAP time series. Finally in Section Hvl we draw our conclusions. 



II. PHASE SHIFTS OF SYNCHRONIZED OSCILLATORS 



A. Winfree model 



The Winfree model is defined by the set of equations {i = 1, ...N) 
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(1) 



It describes a set of 1 coupled non linear oscillators, with coupling constant proportional to k. 9i{t) is the 

phase of the i— th oscillator; {lui} describes a set of natural frequencies taken randomly from a distribution g{Lo). 
P{Oj) is the influence function of the j—th; R{9i) is the sensitivity function giving the response of the i— th oscillator. 
We shall assume below : g{i^) = I/27 for 7 S [wq ~ 77^^o + 7]; 9{^) — otherwise. In the previous equation 
P{0) = 1 + cos6', R{0) ~ — sin6'. The phase diagram of the Winfree model has been recently discussed In 
particular the long-time behavior of the system is characterized by a synchronous dynamics for k and 7 not very 
large. For loq = 1 synchronization occurs for k < 0.77 and 7 < h{K), where the function can be found in Fig. 3 
of Ref. 0; in any case 7 < 0.2. This means that all the oscillators are characterized by a common average frequency 
(or rotation number) pj, = limj^oo Oi{t)/t. The Winfree model can describe different sets of pulse-coupled biological 
osciUators, see e.g. |g. [lOl[lll |. 

We now wish to study the relation between the phase shift 59 of a pair of oscillators and the difference of their 
natural frequencies Slo. We have performed numerical simulations with N = 500 oscillators with different values of 
K and 7 = 0.10 corresponding to the synchronization phase. We have considered times as large a.s t = 1,000. As 
expected there is no dependence on the initial conditions. On the contrary SO is linearly related to Slo as shown in Fig. 
n where we plot pi versus tOi for various values of k (on the left) and S9 versus Sco (on the right). This dependence 
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FIG. 1; On the left: the rotation number p plotted versus uj for 7 — 0.10 and k = 0.35, 0.45, 0.65 (from top to bottom). On 
the right: 56 versus Slo for the same values of 7 and k (larger slopes correspond to smaller values of k). 



can be understood as follows. As A — > 00, the sum over all oscillators in can be replaced by an integral, yielding 

v{e,t,Lo) ^ LO - cr(t) sine (2) 
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where 

p2-rr (•1+7 

a{t) ^ K / {l + cos9)p{9,t,uj)g{uj)dujde . (3) 
Jo Jl-t 

Here p {9, t, w) denotes the density of oscillators with phase 9 at time t. We consider the large t behavior to allow for 
synchronization; moreover we take a temporal averaging over the common period T to get rid of local fluctuations. 
We get 



t+T 

dt a (t) sin 9 (t) (4) 



and consider variations in 



1 



t+T 



= 5lu-- dt a{t) S9{t) cos 9{t) . (5) 
^ Jt 

Since the oscillators are synchronized S9{t) is time-independent for t large enough. Therefore 

k69 /■*+'^ /■ 1'^'^ » » » 

Sui = ~/ dLjg{Lj) / d9{l + cos9)p{9,t,u;) cos9{t) = 

kS9 /■*+'^ kS9 /■*+'^ /■ f'^'^ - 

= / dt cos9{t) + —— / / duj g[uj) / d9p{9,t,u!) cos[9 + 9{t)] + 

+ — dt j dujg{uj) j d9p{9,t,Lj) cos[9-9{t)] . (6) 

The first two terms on the r.h.s. of © vanish since the integrand functions have zero temporal average. We get 
therefore 

Suj^— 69, (7) 
which is the desired linear relation between Suj and 69. The factor A is evaluated as follows: 



I rt+T ^1+7 ^27r 

dwg(cj) / d^p (^^, cos[^- 6i(t)] 
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^"dygf '-21) cos ('-21) . (8) 



Here is probability distribution of 69. It is related to the probability distribution of 6uj hy Both the 6uj 

and the 69 distribution functions are derived from the lo density g{uj). If this density is flat, as assumed here, then 
9 = 9*9, i-e. 

kX / kX^ ^ 



~9m ^ - - (^-1 \69\ . (9) 



In conclusion A is given by solving the equation 

, . 2 ^27 



We notice that in Q there is no dependence on wq; this dependence is in the flrst two terms of JB)) since they vanish 
only in the large N, large t limit. 

In Fig. Elwe report the slope 69/6uj as computed by l|10|l . as a function of k (with 7 = 0.1) on the left and as a 
function of 7 (with k = 0.45) on the right. This curve is independent of w^. We also report results of the numerical 
analysis for two values of ojq These data show a small dependence on 
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FIG. 2; The slope — in the Winfree modeL On the left: — as a function of n for two values of cjq and 7 = 0.1. On the right: 
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50 

— as a function of 7 for two values of ujo and k = 0.45. The curves are independent of ujq. 



B. Kuramoto model 



The analysis of the Kuramoto model produces comparable results. The Kuramoto model is based on the set of 
equations = 1, ...N) 

N 

^,(t)=c^. + ^^sin(0.-0,) . (11) 

The numerical results one gets are similar to those of fig. Q with a linear dependence of 59 on Soj and the rotational 
frequency pi — ojq (we use the same distribution function g(w) as before). The latter results follows from the fact 
that the phases 6i are dynamically pulled toward the the mean phase These results can be compared with an 
analytical treatment by observing that, in this case, instead of © one gets 

6uj=—J^ dt J dujg{uj) J dep{9,t,Lj)cos[9-e{t)] . (12) 

Due to the absence of terms analogous to the first and second terms on the r.h.s. of lO, we expect a better agreement 
between numerical simulations and analytical evaluation. From H12|l we get, instead of 10: 

6uj^KXde. (13) 

with A given by 

These results are reported in fig. [Altogether with the results of the numerical simulations. 



III. SYSTOLIC/DIASTOLIC ARTERIAL PRESSURE RELATION 



A. Phase shifts from arterial pressure data 



Let us consider two time series: xs(t) and xoit), representing systolic and diastolic arterial pressure. Data are 
from a population of 47 normal subjects (mean age+/-SD: 54+/-8 years) who underwent a 10 minutes supine resting 
recording of EGG and noninvasive arterial blood pressure (by the Finapres device), in the laboratory for the assessment 
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FIG. 3: The slope — in the Kuramoto model. On the left: — as a function of k for two values of ljq and 7 — 0.1. On the 
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right: — as a function of 7 for two values of u)o and k = 0.45. 
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of Autonomic Nervous Sytem, S. Maugeri Foundation, Scientific Institute of Montescano, Italy. For each cardiac cycle, 
corresponding values of SAP and DAP were computed and resampled at a frequency of 2 Hz using a cubic spline 
interpolation. In Fig. 2]we report the time series of the systolic arterial pressure xs{t) for one of the subjects examined 
in this study. 
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FIG. 4: The time series xs{t) (on the left) and xoit) (on the right) of the systolic and diastolic arterial pressures for one of 
the subjects examined. 

These data are analyzed by filtering in appropriate frequency bands. We consider here three bands: Very Low 
Frequency (VLF) band: (0.01, 0.04) Hz; Low Frequency (LF) band: (0.04, 0.15) Hz; High Frequency (HF) band: 
(0.15, 0.45) Hz. In a previous paper 4], using Fourier transform methods, occurrence of delays between SAP and 
DAP was investigated, and it was found that DAP anticipates SAP in VLF (delay equal to 2.5 sec )and in LF (0.6 
sees); no significant delay was found in HF. Here we enlarge the statistical population with respect to and evaluate 
the phases of signals by the analytic signal technique, which allows a better estimate. As well known SAP and DAP 
are highly synchronized and our data confirm this expectation. We have used the Hilbert transform method that 
allows to detect phase synchronization in noisy scalar signals [l^. To extract a phase from the signal one considers 
the Hilbert transform of the original time series 



TT J_oo t-T 



(15) 
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where P.V. denotes Cauchy principal value. Then one forms the analytic signal z{t) = x{t) + iy{t) = A(t)e^'^^'^\ 
where A{t) = yjx'^[t) + y^(i) and is the desired phase. To control the possible synchronization of two signals 
xi(t), X2(t) the following procedure is applied: the phases 4'i{t) and 4>2it) are computed and the so called generalized 
phase differences ^n,mit) = [m4>i{t) — 'n4>2{i)]„iod2TT^ with n,m integers, are evaluated. In the present study only 
1 : 1 synchronization has been examined and the two phases <^i(t), 4>2{t) coincide with the phases of the time series 
X]j{t), xs{t). Phase synchronization is characterized by the appearance of peaks in the distribution of the phase 
difference. To evaluate the phase shift we have considered time intervals characterized by a constant phase difference 
between the two series: 

S9 = eoit) ^ es{t) . (16) 

We find 69 > m the VLF band, i.e. in this band diastolic pressure anticipate systolic pressure. Our analysis confirm 
the results of :4] with a different method. On the other hand in the HF band (in basal conditions) the phase shift 
is negative 69 < 0, which means that in this band the systolic pressure signal anticipates the diastolic one. These 
data are reported for all the 47 subjects in fig. Othat shows on the left the VLF band and on the right the HF band 
(we have not reported data in the intermediate region LF band, as they are compatible with 69 — 0). We estimated 
1.76 X 10~^ to be the probability that the phase shifts in the VLF band are sampled from a distribution whose mean is 
less than or equal to zero; analogously 3.0 x 10"^ is the probability that the phase shifts in the HF band are sampled 
from a distribution whose mean is greater than or equal to zero. 




FIG. 5: The phase shifts 56 for all the 47 subjects filtered in the VLF (on the left) and HF (on the right) bands. 

On a physiological basis these results mean that the different sets of oscillators producing the time series have 
different spectral properties. Leaving aside the task of a physiologically-based modelization we now show that the 
results obtained in Section ^ can shed light on this phenomenon. For the present application we use the Winfree 
model. 



B. Interpretation of phase shifts between related oscillatory signals 

We present here a schematic view of phase shifts 69 between the time series xs{t) and xo(t). This picture is only 
qualitative and aims to reproduce the dependence of the sign of 69 on the filter in frequency power spectrum. As 
such, the picture is not realistic and does not offer a physiologically-based model of the time series; nevertheless it can 
shed light on oscillator dynamics underlying them. Let us assume that the two oscillatory signals xs{t) and XD{t) are 
the result of the collective behavior of two sets of oscillators, sets SAP and DAP respectively. We assume that this 
collective behavior produces a Systolic Arterial Pressure (SAP) and Diastolic Arterial Pressure {DAP) time series. 
We assume that the oscillators in the set SAP have natural frequencies in the domain cu e (a — 7, a) 1J(6, 6 -I- 7) , while 
frequencies for the set DAP are in the domain oj e (a, a -I- 7) IJ (6 — 7, 6). We also assume 7 < 1. We will use as 
numerical values a — 1, b = 2 and 7 = 0.1, see FigEl 



On the other hand the two bands A: 1— 7<w<l-|-7, and B: 2 — 7<cj<2 + 7 would model the VLF and HF 
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FIG. 6: The two intervals of natural frequencies A and B. A is on the left and is centered around the value a — 1; B, on the 
right, is centered around fe = 2. We assume that oscillators with frequencies in the band SAP (resp. DAP) produce collectively 
the signal xsit) (resp. XD{t)), see text. 

frequency bands. 

Let us assume that the coupling among the oscillators having natural frequencies in the intervals A and B is 
modelled by the Winfree model, i.e. by eq. (^). However we assume for the coupling 

K^Kij = kH [2j - \uJi - UJjW , (17) 

where H is the Heaviside function. By this choice there is no interaction between oscillators in the two bands, though 
a weak coupling would not alter the qualitative picture. We consider one value of k (k = 0.65 in this case). The 
two sets of oscillators, one centered around the natural frequency luq — 1 (Set A) and the other around luq — 2 (Set 
B) become synchronized around two synchronization frequencies, u>vlf = 0.62 Hz e lovlf = 1-88 Hz see Fig|71(left 
side). Correspondingly, on the right, we have two lines showing a linear dependence between Su) and 59. The two lines 




FIG. 7: On the left: The oscillators of set A, with frequencies centered around a — 1 become synchronized with a frequency 
around ujvlf ~ 0.62 Hz; those of set B (frequencies around 6 = 2) have a synchronization frequency ujhf = 1.88 Hz. On the 
right: The phase shift 59 between any pair of oscillators as a function of the difference Su between the natural frequencies of 
the oscillators in the pairs. The partially overlapping lines refer to the two sets of oscillators A and B, which shows a weak 
dependence of the slope on the natural frequencies. Numerical results refer to A'^ = 1000 oscillators, with k = 0.65. 

are almost completely superimposed, which shows a weak dependence on the average natural frequencies of the two 
sets. The interesting result however is related to our definition of VLF and HF bands. The VLF band is the result 
of the collective behavior of oscillators in the set A. For them ujd > ^^s and therefore, on the basis of the results of 
Section^ 59 = Ooit) — 9s (t) > 0. On the other hand in the HF band, ojs > ^^d and therefore 59 < 0. This simple 
mechanism implies the effect of a change of sign between the two phases when one goes from the very low frequency 
to the high frequency. 

One might wonder if our data also show a casual dependence between the two time series. To address this issue we 
have considered the index S{X\Y) that measures the nonlinear interdependency between two time series X and Y, as 
described in 0. More precisely, from the time series x and y, one reconstructs delay vectors x„ = (a;„, Xn-(m-i)T) 
and y„ = (j/„, yn-(m--i)T)^ where n = 1, N is the time index, m is the embedding dimension, and r denotes the 
time lag. Let rn,j and s„j-, j — 1, A: denote the time indices of the k nearest neighbors of x„ and y„, respectively. 
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For each x„, the mean squared Euchdean distance to its k neighbors is defined as 

i?W(X) = i^(x„-x._f, (18) 

while the Y-conditioned mean squared Euchdean distance is defined by replacing the nearest neighbors by the equal 
time partners of the closest neighbors of y„, 

i?W(X|Y)^i^(x„-x,„j2. (19) 



k 



The interdependence measure is then defined as 



g(X|Y) = 1 y ^" ■ (20) 

A^„tli?i^-)(X|Y) 

S{X\Y) is an asymmetric quantity and the degree of asymmetry is connected to causal relationship between the two 
time series, in other words if S{X\Y) is much greater than S{Y\X) then we may conclude that Y is driving X. On 
the other hand each of these values measures by its size the degree of interdependency, 5 = (resp. 5=1) meaning 
minimal (resp. maximal) interdependency . We evaluated these indexes both on the SAP and DAP time series xs{t), 
Xoit) and on their phases Os{t), Ooit). In both cases the asymmetry was not significant, which means that there is 
no causal relationship between SAP and DAP time series; however the results obtained with the phases are always 
much greater than those obtained with the full signals. Quantitatively, the average values for the HE band are: 
S{xjj\xs) = 4.8 X 10^^, S{xs\xd) = 4.6 x 10~^, which shows a very small asymmetry and, at the same time, a very 
small interdependency. As to the phases, we get S{9D\ds) = 0.899, S{9s\0£)) = 0.901, which on the contrary shows a 
larger interdependency. Similar results are obtained in VLE and LE bands. Besides showing the absence of a causal 
relation, these results confirm that in these systems of oscillators the main source of information on the underlying 
structures resides in the dynamics of the phases. 



IV. CONCLUSIONS 



Our results represent an original analysis of the relation between systolic/diastolic blood pressure, which completes 
previous studies 0| . The measured delays between the oscillatory components of systolic and diastolic blood pressure 
time series, show a change of sign going from low to high frequencies. We have addressed it within the paradigm of 
coupled nonlinear synchronous oscillators. We have shown, using Winfree and Kuramoto models, that once synchro- 
nization is achieved, the phase delay between oscillators is determined by the underlying structure and we have found 
a linear relationship between oscillator phase shifts and the difference of their natural frequencies. We then used 
these results to describe our findings, that confirm that changes in the modulating factors of arterial pressure affect 
differently the systolic and diastolic pressure values 
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